function res = analyze(obj, channelIndex, sampleName, varargin)
%ANALYZE Summary of this function goes here
%   Detailed explanation goes here
    spect = obj.getNodeData(channelIndex, sampleName);
    freq = spect.freq(:).';
    val = spect.value(:).';
    
    p=inputParser;
    p.addParameter('range', minmax(freq), @(x) length(x) == 2 && x(2) > x(1) );
    p.addParameter('func', @semilogy, @(x) ismember(x, {@semilogy, @plot, @loglog}));
    p.addParameter('Npeak', 1, @(x) isscalar(x) && x >= 1);
    p.addParameter('NENBW', 1.5, @(x) isscalar(x) && x > 1.0); % 1.5 for hanning win.
    p.addParameter('uBound', max(val), @(x) isnumeric(x) && isscalar(x) );
    p.addParameter('nBin', 100, @(x) isscalar(x) && x > 1);
    p.addParameter('output', 'all', @(x) ismember(x, {'all', 'peak', 'background'}));
    p.addParameter('factor', 1.0, @(x) isscalar(x) && isnumeric(x) );
    p.addParameter('unit', 'V', @(x) ischar(x));
    p.addParameter('isPlot', true, @islogical);
    p.addParameter('figHdl', [], @(x) isempty(x) || isa(x, 'matlab.ui.Figure'));
    p.parse(varargin{:});
    
    %%
    df = freq(2) - freq(1);
    enbw = p.Results.NENBW*df;
    if ismember(p.Results.output, {'all', 'peak'})       
        rangeIdx = ( freq >= p.Results.range(1) & freq <= p.Results.range(2) );

        [pks, res.locs] = findpeaks(val(rangeIdx), freq(rangeIdx), 'Npeak', p.Results.Npeak, 'SortStr', 'descend');
        res.peaks = pks * sqrt(enbw);
    end
    
    %%
    if ismember(p.Results.output, {'all', 'background'})
        [c, edges]=histcounts(val(val<p.Results.uBound & rangeIdx), p.Results.nBin);
        binSize = edges(2) - edges(1);
        edgesVal = edges(1:end-1)+0.5*binSize;
        [res.background, res.stdBg] = fitGaussian(edgesVal, c);
    end

    %%
    if p.Results.isPlot
        if isempty(p.Results.figHdl)
            figure('Name', 'Spectrum Analyser', 'Position', [100, 100, 1200,800]);
        else
            gcf = p.Results.figHdl;
        end
        ax1 = subplot(1, 4, 1); ax2 = subplot(1, 4, [2 3 4]);
        semilogx(ax1, edgesVal*p.Results.factor, c, '-', [res.background, res.background]*p.Results.factor, minmax(c), 'r--'); grid(ax1, 'on');
        p.Results.func(ax2, freq, val*p.Results.factor, minmax(freq), [res.background, res.background]*p.Results.factor, 'r--'); grid(ax2, 'on'); hold(ax2, 'on');
        
        if p.Results.range(1) > min(freq)
            p.Results.func(ax2, [p.Results.range(1) p.Results.range(1)], ylim(ax2), 'k--');
        end
        if p.Results.range(2) < max(freq)            
            p.Results.func(ax2, [p.Results.range(2) p.Results.range(2)], ylim(ax2), 'k-.');
        end
        
        view(ax1, -90, 90);
        xlim(ax1, ylim(ax2)); xlabel(ax1, ['PSD (' p.Results.unit '\cdot Hz^{-1/2})']); ylabel(ax1, 'Counts');
        
        
        xlim(ax2, minmax(freq));
        yticklabels(ax2, {});
        xlabel(ax2, 'Freuquency (Hz)');
        
        for k = 1:length(res.peaks)
            text(ax2, res.locs(k), res.peaks(k)/sqrt(enbw)*p.Results.factor, sprintf(['\\leftarrow Peak #%d = %3.2e ' p.Results.unit], k,...
                res.peaks(k)*p.Results.factor), 'Color','red','FontSize',16);
        end
        text(ax2, mean(minmax(freq)), (res.background+10.0*res.stdBg)*p.Results.factor, sprintf(['(%3.2e\\pm%2.1e)' p.Results.unit '\\cdot Hz^{-1/2}'],...
            res.background*p.Results.factor, res.stdBg*p.Results.factor), 'Color','red','FontSize',16);
    end
end


function [mu, sigma, amplitude, gof] = fitGaussian(x, y)
    [maxY, maxIdx] = max(y);
    x0 = x(maxIdx);
    [~, idx] = min(abs(y-maxY/exp(1.0)));
    stdGuess = abs(x(max(idx)) - x0);
    
    [xData, yData] = prepareCurveData( x, y );

    % Set up fittype and options.
    ft = fittype( 'gauss1' );  % a1*exp(-((x-b1)/c1)^2)
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower =      [0.9*maxY 0.5*x0 0];
    opts.Upper =      [1.1*maxY 2.0*x0 Inf];
    opts.StartPoint = [    maxY     x0 stdGuess];

    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );
    mu = fitresult.b1;
    sigma = fitresult.c1/sqrt(2.0);
    amplitude = fitresult.a1;
end